library(magrittr)
library(stringr)
library(plyr)
library(tidyverse)


rm(list=ls())
home = 'C:/Users/Jason/Dropbox/VNA_Responsiveness/Analysis/JOP-dataverse/'


permutation.results = paste0(home, 'RI-analyses.Rds') %>%
  readRDS
experimental.results = paste0(home, 'experimental-analyses.Rds') %>%
  readRDS


Random = subset(permutation.results, 
                Result %in% c('tab2col1','tab2col2','tab2col6','placebo') & 
                  !term %in% c('FullTime','CentNom','Competitive'))
Actual = ddply(.data=Random, 
               .variables=c('Result','term'), 
               .fun=function(Coefficient) {
                 X = Coefficient$estimate
                 qest = ecdf(X)
                 Y = subset(experimental.results, 
                            Result==Coefficient$Result[1] & term==Coefficient$term[1])$estimate
                 Q = mapvalues(x=as.integer(round(1000*qest(Y))), 
                               from=seq(0L, 1000L, 1L),
                               to=paste0('(', 
                                         str_pad(string=c(seq(from=0L, 
                                                              to=499L, 
                                                              by=1L)/500L, 
                                                          0.500, 
                                                          seq(from=499L, 
                                                              to=0L, 
                                                              by=-1L)/500L), 
                                                 width=5, 
                                                 side='right', 
                                                 pad=0), 
                                         ')'),
                               warn_missing=F)
                 data.frame(estimate=Y, 
                            Significance=Q, 
                            stringsAsFactors=F)
               }, .inform=T) %>%
  mutate(Term=mapvalues(x=term, 
                        from=c('Citizen',
                               'Citizen:Education',
                               'Firm',
                               'Firm:Education',
                               'Education',
                               '(Intercept)'),
                        to=6:1,
                        warn_missing=F),
         Term=as.integer(Term),
         Result=mapvalues(x=Result, 
                          from=c('tab2col1',
                                 'tab2col2',
                                 'tab2col6',
                                 'placebo'), 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F),
         Result=factor(x=Result,
                       levels=c('1) All delegates\nn = 293',
                                '2) All delegates + covariates\nn = 293',
                                '3) No central nominees\nn = 250',
                                '4) Placebo test\nn = 879')))
Random = Random %>%
  mutate(Term=mapvalues(x=term, 
                        from=c('Citizen',
                               'Citizen:Education',
                               'Firm',
                               'Firm:Education',
                               'Education',
                               '(Intercept)'),
                        to=6:1,
                        warn_missing=F),
         Term=as.integer(Term),
         Result=mapvalues(x=Result, 
                          from=c('tab2col1',
                                 'tab2col2',
                                 'tab2col6',
                                 'placebo'), 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F),
         Result=factor(x=Result,
                       levels=c('1) All delegates\nn = 293',
                                '2) All delegates + covariates\nn = 293',
                                '3) No central nominees\nn = 250',
                                '4) Placebo test\nn = 879')))
Bounds = (max(abs(c(Actual$estimate, Random$estimate)))+0.1)
xLabels = seq(from=(-1)*max(floor(Bounds), 1), 
              to=max(floor(Bounds), 1), 
              by=1L)
Effects = paste(str_pad(string=str_extract(string=Actual$estimate, 
                                           pattern='.*(?=\\.)'), 
                        width=2, 
                        side='left', 
                        pad=' '),
                str_pad(string=str_extract(string=round(Actual$estimate, 3), 
                                           pattern='(?<=\\.).*'), 
                        width=3, 
                        side='right', 
                        pad='0'), 
                sep='.') %>%
  as.data.frame(stringsAsFactors=F) %>%
  set_colnames('Label') %>%
  mutate(estimate=xLabels[1]+0.05, 
         Term=Actual$Term, 
         Result=Actual$Result,
         Significance=Actual$Significance,
         Label=str_c(Label, '\n', Significance))
tmp = subset(Random, Term==6L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Citizen = ggplot_build(tmp)$data[[1]] %>% 
  mutate(scaled=(scaled*0.8)+6, ymin=ymin+6,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
tmp = subset(Random, Term==5L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Citizen.Education = ggplot_build(tmp)$data[[1]] %>% 
  mutate(scaled=(scaled*0.8)+5, ymin=ymin+5,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
tmp = subset(Random, Term==4L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Firm = ggplot_build(tmp)$data[[1]] %>% 
  mutate(scaled=(scaled*0.8)+4, ymin=ymin+4,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
tmp = subset(Random, Term==3L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Firm.Education = ggplot_build(tmp)$data[[1]] %>%
  mutate(scaled=(scaled*0.8)+3, ymin=ymin+3,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
tmp = subset(Random, Term==2L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Education = ggplot_build(tmp)$data[[1]] %>%
  mutate(scaled=(scaled*0.8)+2, ymin=ymin+2,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
tmp = subset(Random, Term==1L) %>%
  ggplot(aes(x=estimate)) +
  geom_density() +
  facet_wrap(~Result, nrow=2, drop=F)
Constant = ggplot_build(tmp)$data[[1]] %>%
  mutate(scaled=(scaled*0.8)+1, ymin=ymin+1,
         Result=mapvalues(x=PANEL, 
                          from=1:4, 
                          to=c('1) All delegates\nn = 293',
                               '2) All delegates + covariates\nn = 293',
                               '3) No central nominees\nn = 250',
                               '4) Placebo test\nn = 879'), 
                          warn_missing=F))
g = ggplot(Random, aes(x=estimate, y=Term)) +
  geom_vline(xintercept=0, linetype=2, color='red') +
  geom_point(size=1, alpha=0.01) +
  geom_text(data=Effects, 
            aes(x=estimate, y=Term, label=Label), 
            size=2.5, alpha=1, color='black', hjust=0, inherit.aes=F) +
  geom_ribbon(data=Citizen, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_ribbon(data=Citizen.Education, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_ribbon(data=Firm, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_ribbon(data=Firm.Education, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_ribbon(data=Education, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_ribbon(data=Constant, 
              aes(x=x, ymin=ymin, ymax=scaled), 
              fill='black', color='black', alpha=0.25, size=0.25, inherit.aes=F) +
  geom_segment(data=Actual, 
               aes(x=estimate, xend=estimate, y=Term, yend=Term+0.8), 
               size=0.25, alpha=1, linetype=1, inherit.aes=F) +
  geom_point(data=Actual, 
             aes(x=estimate, y=Term), 
             fill='white', size=2, alpha=1, shape=21, inherit.aes=F) +
  scale_x_continuous(breaks=xLabels) +
  scale_y_continuous(breaks=6:1, 
                     labels=c('Citizen',
                              'Citizen \u00D7 EducationLaw',
                              'Firm',
                              'Firm \u00D7 EducationLaw',
                              'EducationLaw',
                              'Constant')) +
  facet_wrap(~Result, nrow=1) +
  labs(x=NULL, y=NULL) +
  guides(fill=guide_legend(override.aes=list(size=3))) +
  coord_cartesian(xlim=range(xLabels), ylim=c(0.6, 7), expand=F) +
  theme_minimal() +
  theme(axis.text=element_text(size=8, color='black'),
        strip.text=element_text(size=9),
        legend.title=element_text(size=9),
        legend.text=element_text(size=8),
        legend.position='bottom',
        plot.margin=unit(c(0.1, 0.2, 0.1, 0.1), 'in'),
        panel.spacing.x=unit(0.2, 'in'),
        panel.border=element_rect(color='black', fill=NA))   


ggsave('figure-05.png', g, path=home, width=9, height=3, units='in')
ggsave('figure-05.tiff', g, path=home, width=9, height=3, units='in')
ggsave('figure-05.eps', g, path=home, width=9, height=3, units='in', device=cairo_ps)
